home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / rng / test.c < prev    next >
Encoding:
C/C++ Source or Header  |  2002-04-18  |  13.6 KB  |  516 lines
  1. /* rng/test.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 James Theiler, Brian Gough
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. #include <config.h>
  21. #include <stdlib.h>
  22. #include <stdio.h>
  23. #include <math.h>
  24. #include <gsl/gsl_rng.h>
  25. #include <gsl/gsl_test.h>
  26. #include <gsl/gsl_ieee_utils.h>
  27.  
  28. void rng_test (const gsl_rng_type * T, unsigned long int seed, unsigned int n,
  29.            unsigned long int result);
  30. void rng_float_test (const gsl_rng_type * T);
  31. void generic_rng_test (const gsl_rng_type * T);
  32. void rng_state_test (const gsl_rng_type * T);
  33. void rng_parallel_state_test (const gsl_rng_type * T);
  34. int rng_max_test (gsl_rng * r, unsigned long int *kmax, unsigned long int ran_max) ;
  35. int rng_min_test (gsl_rng * r, unsigned long int *kmin, unsigned long int ran_min, unsigned long int ran_max) ;
  36. int rng_sum_test (gsl_rng * r, double *sigma);
  37. int rng_bin_test (gsl_rng * r, double *sigma);
  38.  
  39. #define N  10000
  40. #define N2 200000
  41.  
  42. int
  43. main (void)
  44. {
  45.   const gsl_rng_type ** rngs = gsl_rng_types_setup();  /* get all rng types */
  46.  
  47.   const gsl_rng_type ** r ;
  48.  
  49.   gsl_ieee_env_setup ();
  50.  
  51.   gsl_rng_env_setup ();
  52.  
  53.   /* specific tests of known results for 10000 iterations with seed = 1 */
  54.  
  55.   rng_test (gsl_rng_rand, 1, 10000, 1910041713);
  56.   rng_test (gsl_rng_randu, 1, 10000, 1623524161);
  57.   rng_test (gsl_rng_cmrg, 1, 10000, 719452880);
  58.   rng_test (gsl_rng_minstd, 1, 10000, 1043618065);
  59.   rng_test (gsl_rng_mrg, 1, 10000, 2064828650);
  60.   rng_test (gsl_rng_taus, 1, 10000, 2733957125UL);
  61.   rng_test (gsl_rng_transputer, 1, 10000, 1244127297UL);
  62.   rng_test (gsl_rng_vax, 1, 10000, 3051034865UL);
  63.  
  64.   rng_test (gsl_rng_borosh13, 1, 10000, 2513433025UL);
  65.   rng_test (gsl_rng_fishman18, 1, 10000, 1626638977UL);
  66.   rng_test (gsl_rng_fishman2x, 1, 10000, 1158622502UL);
  67.   rng_test (gsl_rng_knuthran2, 1, 10000, 1182784902UL);
  68.   rng_test (gsl_rng_knuthran, 310952, 1009 * 2009 + 1, 461390032);
  69.   rng_test (gsl_rng_lecuyer21, 1, 10000, 2006618587UL);
  70.   rng_test (gsl_rng_waterman14, 1, 10000, 3776680385UL);
  71.  
  72.   /* specific tests of known results for 10000 iterations with seed = 6 */
  73.  
  74.   rng_test (gsl_rng_coveyou, 6, 10000, 1416754246UL);
  75.   rng_test (gsl_rng_fishman20, 6, 10000, 1811577350UL);
  76.  
  77.   /* FIXME: the ranlux tests below were made by running the fortran code and
  78.      getting the expected value from that. An analytic calculation
  79.      would be preferable. */
  80.  
  81.   rng_test (gsl_rng_ranlux, 314159265, 10000, 12077992);
  82.   rng_test (gsl_rng_ranlux389, 314159265, 10000, 165942);
  83.  
  84.   rng_test (gsl_rng_ranlxs0, 1, 10000, 11904320);
  85.   /* 0.709552764892578125 * ldexp(1.0,24) */
  86.  
  87.   rng_test (gsl_rng_ranlxs1, 1, 10000, 8734328);
  88.   /* 0.520606517791748047 * ldexp(1.0,24) */
  89.  
  90.   rng_test (gsl_rng_ranlxs2, 1, 10000, 6843140); 
  91.   /* 0.407882928848266602 * ldexp(1.0,24) */
  92.  
  93.   rng_test (gsl_rng_ranlxd1, 1, 10000, 1998227290UL);
  94.   /* 0.465248546261094020 * ldexp(1.0,32) */
  95.  
  96.   rng_test (gsl_rng_ranlxd2, 1, 10000, 3949287736UL);
  97.   /* 0.919515205581550532 * ldexp(1.0,32) */
  98.  
  99.   /* FIXME: the tests below were made by running the original code in
  100.      the ../random directory and getting the expected value from
  101.      that. An analytic calculation would be preferable. */
  102.  
  103.   rng_test (gsl_rng_slatec, 1, 10000, 45776);
  104.   rng_test (gsl_rng_uni, 1, 10000, 9214);
  105.   rng_test (gsl_rng_uni32, 1, 10000, 1155229825);
  106.   rng_test (gsl_rng_zuf, 1, 10000, 3970);
  107.  
  108.   /* The tests below were made by running the original code and
  109.      getting the expected value from that. An analytic calculation
  110.      would be preferable. */
  111.  
  112.   rng_test (gsl_rng_r250, 1, 10000, 1100653588);
  113.   rng_test (gsl_rng_mt19937, 4357, 1000, 1030650439);
  114.   rng_test (gsl_rng_tt800, 0, 10000, 2856609219UL);
  115.  
  116.   rng_test (gsl_rng_ran0, 0, 10000, 1115320064);
  117.   rng_test (gsl_rng_ran1, 0, 10000, 1491066076);
  118.   rng_test (gsl_rng_ran2, 0, 10000, 1701364455);
  119.   rng_test (gsl_rng_ran3, 0, 10000, 186340785);
  120.  
  121.   rng_test (gsl_rng_ranmar, 1, 10000, 14428370);
  122.  
  123.   rng_test (gsl_rng_rand48, 0, 10000, 0xDE095043UL);
  124.   rng_test (gsl_rng_rand48, 1, 10000, 0xEDA54977UL);
  125.  
  126.   rng_test (gsl_rng_random_glibc2, 0, 10000, 1908609430);
  127.   rng_test (gsl_rng_random8_glibc2, 0, 10000, 1910041713);
  128.   rng_test (gsl_rng_random32_glibc2, 0, 10000, 1587395585);
  129.   rng_test (gsl_rng_random64_glibc2, 0, 10000, 52848624);
  130.   rng_test (gsl_rng_random128_glibc2, 0, 10000, 1908609430);
  131.   rng_test (gsl_rng_random256_glibc2, 0, 10000, 179943260);
  132.  
  133.   rng_test (gsl_rng_random_bsd, 0, 10000, 1457025928);
  134.   rng_test (gsl_rng_random8_bsd, 0, 10000, 1910041713);
  135.   rng_test (gsl_rng_random32_bsd, 0, 10000, 1663114331);
  136.   rng_test (gsl_rng_random64_bsd, 0, 10000, 864469165);
  137.   rng_test (gsl_rng_random128_bsd, 0, 10000, 1457025928);
  138.   rng_test (gsl_rng_random256_bsd, 0, 10000, 1216357476);
  139.  
  140.   rng_test (gsl_rng_random_libc5, 0, 10000, 428084942);
  141.   rng_test (gsl_rng_random8_libc5, 0, 10000, 1910041713);
  142.   rng_test (gsl_rng_random32_libc5, 0, 10000, 1967452027);
  143.   rng_test (gsl_rng_random64_libc5, 0, 10000, 2106639801);
  144.   rng_test (gsl_rng_random128_libc5, 0, 10000, 428084942);
  145.   rng_test (gsl_rng_random256_libc5, 0, 10000, 116367984);
  146.  
  147.   rng_test (gsl_rng_ranf, 0, 10000, 2152890433UL);
  148.   rng_test (gsl_rng_ranf, 2, 10000, 339327233);
  149.  
  150.   /* Test constant relationship between int and double functions */
  151.  
  152.   for (r = rngs ; *r != 0; r++)
  153.     rng_float_test (*r);
  154.  
  155.   /* Test save/restore functions */
  156.  
  157.   for (r = rngs ; *r != 0; r++)
  158.     rng_state_test (*r);
  159.  
  160.   for (r = rngs ; *r != 0; r++)
  161.     rng_parallel_state_test (*r);
  162.  
  163.   /* generic statistical tests (these are just to make sure that we
  164.      don't get any crazy results back from the generator, i.e. they
  165.      aren't a test of the algorithm, just the implementation) */
  166.  
  167.   for (r = rngs ; *r != 0; r++)
  168.     generic_rng_test (*r);
  169.  
  170.   exit (gsl_test_summary ());
  171. }
  172.  
  173.  
  174. void
  175. rng_test (const gsl_rng_type * T, unsigned long int seed, unsigned int n,
  176.       unsigned long int result)
  177. {
  178.   gsl_rng *r = gsl_rng_alloc (T);
  179.   unsigned int i;
  180.   unsigned long int k = 0;
  181.   int status;
  182.  
  183.   if (seed != 0)
  184.     {
  185.       gsl_rng_set (r, seed);
  186.     }
  187.  
  188.   for (i = 0; i < n; i++)
  189.     {
  190.       k = gsl_rng_get (r);
  191.     }
  192.  
  193.   status = (k != result);
  194.   gsl_test (status, "%s, %u steps (%u observed vs %u expected)",
  195.         gsl_rng_name (r), n, k, result);
  196.  
  197.   gsl_rng_free (r);
  198. }
  199.  
  200. void
  201. rng_float_test (const gsl_rng_type * T)
  202. {
  203.   gsl_rng *ri = gsl_rng_alloc (T);
  204.   gsl_rng *rf = gsl_rng_alloc (T);
  205.  
  206.   double u, c ; 
  207.   unsigned int i;
  208.   unsigned long int k = 0;
  209.   int status = 0 ;
  210.  
  211.   do 
  212.     {
  213.       k = gsl_rng_get (ri);
  214.       u = gsl_rng_get (rf);
  215.     } 
  216.   while (k == 0) ;
  217.  
  218.   c = k / u ;
  219.   for (i = 0; i < N2; i++)
  220.     {
  221.       k = gsl_rng_get (ri);
  222.       u = gsl_rng_get (rf);
  223.       if (c*k != u)
  224.     {
  225.       status = 1 ;
  226.       break ;
  227.     }
  228.     }
  229.  
  230.   gsl_test (status, "%s, ratio of int to double (%g observed vs %g expected)",
  231.         gsl_rng_name (ri), c, k/u);
  232.  
  233.   gsl_rng_free (ri);
  234.   gsl_rng_free (rf);
  235. }
  236.  
  237.  
  238. void
  239. rng_state_test (const gsl_rng_type * T)
  240. {
  241.   unsigned long int test_a[N], test_b[N];
  242.  
  243.   int i;
  244.  
  245.   gsl_rng *r = gsl_rng_alloc (T);
  246.   gsl_rng *r_save = gsl_rng_alloc (T);
  247.  
  248.   for (i = 0; i < N; ++i)
  249.     {
  250.       gsl_rng_get (r);    /* throw away N iterations */
  251.     }
  252.  
  253.   gsl_rng_memcpy (r_save, r);    /* save the intermediate state */
  254.  
  255.   for (i = 0; i < N; ++i)
  256.     {
  257.       test_a[i] = gsl_rng_get (r);
  258.     }
  259.  
  260.   gsl_rng_memcpy (r, r_save);    /* restore the intermediate state */
  261.   gsl_rng_free (r_save);
  262.  
  263.   for (i = 0; i < N; ++i)
  264.     {
  265.       test_b[i] = gsl_rng_get (r);
  266.     }
  267.  
  268.   {
  269.     int status = 0;
  270.     for (i = 0; i < N; ++i)
  271.       {
  272.     status |= (test_b[i] != test_a[i]);
  273.       }
  274.     gsl_test (status, "%s, random number state consistency",
  275.           gsl_rng_name (r));
  276.   }
  277.  
  278.   gsl_rng_free (r);
  279. }
  280.  
  281.  
  282. void
  283. rng_parallel_state_test (const gsl_rng_type * T)
  284. {
  285.   unsigned long int test_a[N], test_b[N];
  286.   unsigned long int test_c[N], test_d[N];
  287.   double test_e[N], test_f[N];
  288.  
  289.   int i;
  290.  
  291.   gsl_rng *r1 = gsl_rng_alloc (T);
  292.   gsl_rng *r2 = gsl_rng_alloc (T);
  293.  
  294.   for (i = 0; i < N; ++i)
  295.     {
  296.       gsl_rng_get (r1);        /* throw away N iterations */
  297.     }
  298.  
  299.   gsl_rng_memcpy (r2, r1);        /* save the intermediate state */
  300.  
  301.   for (i = 0; i < N; ++i)
  302.     {
  303.       /* check that there is no hidden state intermixed between r1 and r2 */
  304.       test_a[i] = gsl_rng_get (r1);    
  305.       test_b[i] = gsl_rng_get (r2);
  306.       test_c[i] = gsl_rng_uniform_int (r1, 1234);    
  307.       test_d[i] = gsl_rng_uniform_int (r2, 1234);
  308.       test_e[i] = gsl_rng_uniform (r1);    
  309.       test_f[i] = gsl_rng_uniform (r2);
  310.     }
  311.  
  312.   {
  313.     int status = 0;
  314.     for (i = 0; i < N; ++i)
  315.       {
  316.     status |= (test_b[i] != test_a[i]);
  317.     status |= (test_c[i] != test_d[i]);
  318.     status |= (test_e[i] != test_f[i]);
  319.       }
  320.     gsl_test (status, "%s, parallel random number state consistency",
  321.           gsl_rng_name (r1));
  322.   }
  323.  
  324.   gsl_rng_free (r1);
  325.   gsl_rng_free (r2);
  326.  
  327. }
  328.  
  329. void
  330. generic_rng_test (const gsl_rng_type * T)
  331. {
  332.   gsl_rng *r = gsl_rng_alloc (T);
  333.   const char *name = gsl_rng_name (r);
  334.   unsigned long int kmax = 0, kmin = 1000;
  335.   double sigma = 0;
  336.   const unsigned long int ran_max = gsl_rng_max (r);
  337.   const unsigned long int ran_min = gsl_rng_min (r);
  338.  
  339.   int status = rng_max_test (r, &kmax, ran_max);
  340.  
  341.   gsl_test (status,
  342.         "%s, observed vs theoretical maximum (%lu vs %lu)",
  343.         name, kmax, ran_max);
  344.  
  345.   status = rng_min_test (r, &kmin, ran_min, ran_max);
  346.  
  347.   gsl_test (status,
  348.         "%s, observed vs theoretical minimum (%lu vs %lu)",
  349.         name, kmin, ran_min);
  350.  
  351.   status = rng_sum_test (r, &sigma);
  352.  
  353.   gsl_test (status,
  354.         "%s, sum test within acceptable sigma (observed %.2g sigma)",
  355.         name, sigma);
  356.  
  357.   status = rng_bin_test (r, &sigma);
  358.  
  359.   gsl_test (status,
  360.         "%s, bin test within acceptable chisq (observed %.2g sigma)",
  361.         name, sigma);
  362.  
  363.   gsl_rng_set (r, 1);    /* set seed to 1 */
  364.   status = rng_max_test (r, &kmax, ran_max);
  365.  
  366.   gsl_rng_set (r, 1);    /* set seed to 1 */
  367.   status |= rng_min_test (r, &kmin, ran_min, ran_max);
  368.  
  369.   gsl_rng_set (r, 1);    /* set seed to 1 */
  370.   status |= rng_sum_test (r, &sigma);
  371.  
  372.   gsl_rng_set (r, 12345);    /* set seed to a "typical" value */
  373.   status |= rng_max_test (r, &kmax, ran_max);
  374.  
  375.   gsl_rng_set (r, 12345);    /* set seed to a "typical" value */
  376.   status |= rng_min_test (r, &kmin, ran_min, ran_max);
  377.  
  378.   gsl_rng_set (r, 12345);    /* set seed to a "typical" value */
  379.   status |= rng_sum_test (r, &sigma);
  380.  
  381.   gsl_test (status, "%s, maximum and sum tests for non-default seeds", name);
  382.  
  383.   gsl_rng_free (r);
  384. }
  385.  
  386. int
  387. rng_max_test (gsl_rng * r, unsigned long int *kmax, unsigned long int ran_max)
  388. {
  389.   unsigned long int actual_uncovered;
  390.   double expect_uncovered;
  391.   int status;
  392.   unsigned long int max = 0;
  393.   int i;
  394.  
  395.   for (i = 0; i < N2; ++i)
  396.     {
  397.       unsigned long int k = gsl_rng_get (r);
  398.       if (k > max)
  399.     max = k;
  400.     }
  401.  
  402.   *kmax = max;
  403.  
  404.   actual_uncovered = ran_max - max;
  405.   expect_uncovered = (double) ran_max / (double) N2;
  406.  
  407.   status = (max > ran_max) || (actual_uncovered > 7 * expect_uncovered) ;
  408.  
  409.   return status;
  410. }
  411.  
  412. int
  413. rng_min_test (gsl_rng * r, unsigned long int *kmin, 
  414.           unsigned long int ran_min, unsigned long int ran_max)
  415. {
  416.   unsigned long int actual_uncovered;
  417.   double expect_uncovered;
  418.   int status;
  419.   unsigned long int min = 1000000000UL;
  420.   int i;
  421.  
  422.   for (i = 0; i < N2; ++i)
  423.     {
  424.       unsigned long int k = gsl_rng_get (r);
  425.       if (k < min)
  426.     min = k;
  427.     }
  428.  
  429.   *kmin = min;
  430.  
  431.   actual_uncovered = min - ran_min;
  432.   expect_uncovered = (double) ran_max / (double) N2;
  433.  
  434.   status = (min < ran_min) || (actual_uncovered > 7 * expect_uncovered);
  435.  
  436.   return status;
  437. }
  438.  
  439. int
  440. rng_sum_test (gsl_rng * r, double *sigma)
  441. {
  442.   double sum = 0;
  443.   int i, status;
  444.  
  445.   for (i = 0; i < N2; ++i)
  446.     {
  447.       double x = gsl_rng_uniform (r) - 0.5;
  448.       sum += x;
  449.     }
  450.  
  451.   sum /= N2;
  452.  
  453.   /* expect the average to have a variance of 1/(12 n) */
  454.  
  455.   *sigma = sum * sqrt (12.0 * N2);
  456.  
  457.   /* more than 3 sigma is an error */
  458.   
  459.   status = (fabs (*sigma) > 3 || fabs(*sigma) < 0.003);
  460.  
  461.   if (status) {
  462.       fprintf(stderr,"sum=%g, sigma=%g\n",sum,*sigma);
  463.   }
  464.  
  465.   return status;
  466. }
  467.  
  468. #define BINS 17
  469. #define EXTRA 10
  470. int
  471. rng_bin_test (gsl_rng * r, double *sigma)
  472. {
  473.   int count[BINS+EXTRA];
  474.   double chisq = 0;
  475.   int i, status;
  476.  
  477.   for (i = 0; i < BINS+EXTRA; i++)
  478.       count[i] = 0 ;
  479.  
  480.  
  481.   for (i = 0; i < N2; i++)
  482.     {
  483.       int j = gsl_rng_uniform_int (r, BINS);
  484.       count[j]++ ;
  485.     }
  486.  
  487.   chisq = 0 ;
  488.   for (i = 0; i < BINS; i++)
  489.     {
  490.       double x = (double)N2/(double)BINS ;
  491.       double d = (count[i] - x) ;
  492.       chisq += (d*d) / x;
  493.     }
  494.  
  495.   *sigma = sqrt(chisq/BINS) ;
  496.  
  497.   /* more than 3 sigma is an error */
  498.   
  499.   status = (fabs (*sigma) > 3 || fabs(*sigma) < 0.003);
  500.  
  501.   for (i = BINS; i < BINS+EXTRA; i++)
  502.     {
  503.       if (count[i] != 0)
  504.     {
  505.       status = 1 ;
  506.       gsl_test (status, 
  507.             "%s, wrote outside range in bin test "
  508.             "(%d observed vs %d expected)",
  509.             gsl_rng_name(r), i, BINS - 1);
  510.     }
  511.     }
  512.  
  513.   return status;
  514. }
  515.  
  516.